The Inverse Real Discrete Fourier Transform¶

In [1]:
%load_ext autoreload
%autoreload 2
import numpy as np
import matplotlib.pyplot as plt
import IPython.display as ipd

If we have an "analysis signal" $x[n]$ with $N$ samples, we define the Discrete Fourier Transform (DFT) by taking dot products with "probe signals," which are sines and cosines that go through an integer number of periods over the interval of $N$ samples.

Forward Discrete Fourier Transform¶

$c[k] = \sum_{n=0}^{N-1} x[n] \cos(2 \pi k n / N) / (N/2)$¶

$s[k] = \sum_{n=0}^{N-1} x[n] \sin(2 \pi k n / N) / (N/2)$¶

In [2]:
def get_dft(x):
    N = len(x)
    t = np.linspace(0, 1, N+1)[0:N]
    n_freqs = int(np.ceil(N/2))
    cos_sums = np.zeros(n_freqs)
    sin_sums = np.zeros(n_freqs)
    for i in range(n_freqs):
        c = np.cos(2*np.pi*i*t)
        s = np.sin(2*np.pi*i*t)
        cos_sums[i] = np.sum(c*x)/(N/2)
        sin_sums[i] = np.sum(s*x)/(N/2)
    return cos_sums, sin_sums
In [3]:
N = 100
n = np.arange(N)
x = 2*np.cos(2*np.pi*(3)*n/N)
x += np.sin(2*np.pi*(5)*(n-2)/N)
c, s = get_dft(x)

plt.figure(figsize=(12, 3))
plt.plot(x)
plt.title("Original Signal")

plt.figure(figsize=(8, 4))
plt.subplot(121)
plt.stem(c[0:10])
plt.title("Cosine Components")
plt.subplot(122)
plt.stem(s[0:10])
plt.title("Sine Components")
Out[3]:
Text(0.5, 1.0, 'Sine Components')

Inverse Discrete Fourier Transform¶

The inverse Discrete Fourier Transform is going backwards; that is, given the result of doing dot products with all cosine and sine probes, go back and reconstruct the analysis signal. This can be done with the following formula and code

$ x[n] = \sum_{k = 0}^{N/2} c[k] \cos(2 \pi k n / N) + s[k] \sin(2 \pi k n / N)$¶

In [4]:
def get_idft(c, s):
    """
    Given the cosine and sine components, reconstruct the signal x
    in the time domain
    """
    N = 2*len(c)
    n = np.arange(N)
    x = np.zeros(N)
    n_freqs = len(c)
    for k in range(n_freqs):
        x += c[k]*np.cos(2*np.pi*(k/N)*n)
        x += s[k]*np.sin(2*np.pi*(k/N)*n)
    return x

We then check to see if our answer agrees with the signal we started with, and indeed it does!

In [5]:
x_recon = get_idft(c, s)
plt.figure(figsize=(12, 4))
plt.plot(x_recon)
plt.plot(x, linestyle='--')
plt.legend(["Reconstructed Signal", "Original Signal"])
Out[5]:
<matplotlib.legend.Legend at 0x7f9af6111130>

Fast(er) Risset Beats¶

In [6]:
from risset import *

def do_risset_faster(filename, tune_length, freqs_per_note, sr):
    """
    Implement the faster version of Risset beats that aggregates
    duplicate frequencies into a sine and cosine term, and then
    uses the fast inverse Discrete Fourier Transform to reconstruct
    the time signal

    Parameters
    ----------
    filename: string
        Path to file with the tune
    tune_length: float
        Length, in seconds, of the tune
    freqs_per_note: int
        Number of frequencies to use for each note
    sr: int
        The sample rate of the entire piece
    """
    ps, times = load_tune(filename, tune_length)
    N = int(sr*tune_length)
    c = np.zeros(N) # Analagous to dictionaries
    s = np.zeros(N) # but using units of cycles per interval (instead of hz, cycles/sec)
    for p, shift in zip(ps, times):
        fhz = 440*(2**(p/12)) # In cycles/second (hz)
        k_center = int(fhz*N/sr)
        # Goal with risset beats to create a beat frequency that happens once
        # over our tune interval
        # Frequency k_center + 1 cycles/interval, difference is 1 cycle/interval
        for dk in range(-freqs_per_note//2, freqs_per_note//2+1):
            k = k_center + dk
            f = k*sr/N
            c[k] += np.cos(2*np.pi*f*shift)
            s[k] += np.sin(2*np.pi*f*shift)
    
    ## Setting up format to use the "Fast Fourier Transform," which is an O(N log N) algorithm
    ## whereas the "brute force" version we implemented is O(N^2)
    M = N//2
    if N%2 == 1:
        M = (N-1)//2
    else:
        M = (N-2)//2
    c[-M::] = c[1:M+1][::-1]
    s[-M::] = -s[1:M+1][::-1]
    return np.real(np.fft.ifft(c - 1j*s))


import time
tune_length = 70
freqs_per_note = 500
sr = 44100
tic = time.time()
x = do_risset_faster("Tunes/wanna.txt", tune_length, freqs_per_note, sr)
print("Elapsed Time: {:.3f}".format(time.time()-tic))
ipd.Audio(x, rate=sr)
Elapsed Time: 0.241
Out[6]:
Your browser does not support the audio element.

Let's look at the difference

In [7]:
N = 44100*70
# N^2 / N log2(N),   N/log2(N)
N/np.log2(N)
Out[7]:
143196.6024188806
In [ ]: